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This  report  describes  a  method  for  computing 
solutions  to  partial  differential  equations  in 
unbounded  regions.  The  differential  equations  and 
unbounded  regions  arise  in  problems  such  as  under¬ 
water  acoustic  and  magnetic  or  electric  fields  sur¬ 
rounding  vehicles  or  instruments  in  the  sea  or  air. 

The  method  couples  finite  element  and  boundary  in¬ 
tegral  equation  solutions  to  compute  shock  response 
of  submarine  hulls,  acoustic  scattering  due  to  elas¬ 
tic  structures,  electromagnetic  signatures  of  ships 
and  submarines,  etc.  The  finite  element  method  is 
used  in  the  bounded  region  occupied  by  the  vehicle 
or  instrument.  In  this  region  there  may  be  acoustic, 
magnetic  or  electric  sources  and  the  material  prop¬ 
erties  may  be  nonuniform.  The  boundary  integral 
equation  method  is  used  in  the  unbounded  region  occu¬ 
pied  by  the  surrounding  sea  or  atmosphere.  This 
method  requires  that  the  unbounded  region  be  assumed 
to  be  free  of  sources  and  have  uniform  material  prop¬ 
erties.  The  coupled  method  produced  good  results  for 
a  sample  problem  in  which  the  temperature  distribution 
was  computed  for  a  two-dimensional  steady  state  heat 
flow.  / 

r-. 

ADMINISTRATIVE  INFORMATION 

This  work  was  funded  by  the  Naval  Sea  Systems  Coanand  (S03R22)  under  Task 
Area  SR0140301,  Task  15321  and  Work  Unit  1808-010. 

I.  INTRODUCTION 

Several  problems  of  interest  to  the  Navy  Involve  a  vehicle  surrounded  by 
a  medium  of  large  extent.  Two  examples  are: 

a  Sound  radiation  or  scattering  from  a  submerged  structure.  The 
structure  may  have  a  complicated  shape,  including  appendages, 
and  the  surrounding  water  is  assumed  to  have  uniform  acoustical 
properties . 

•  Electric  or  magnetic  fields  surrounding  a  ship.  The  ship  may 
have  an  Irregular  shape  and  may  have  machinery  and  materials  with 
differing  magnetic  properties  and  machinery  that  influences  the 
fields,  but  the  surrounding  air  or  water  is  assumed  to  have  uni¬ 
form  electric  and  magnetic  properties. 


The  characteristic  feature  of  such  problems  is  that  the  vehicle  is  relatively 
small  and  may  have  a  complex  configuration,  while  the  surrounding  medium  has  a 
large  extent  and  is  assumed  to  have  uniform  properties.  These  problems  present 
two  contrasting  requirements.  There  is  a  small  inner  region  that  requires  a 
detailed  representation  and  an  outer  region  that  does  not  need  detailed  repre¬ 
sentation  but  has  a  large  extent.  The  outer  region  is  often  so  large  that  its 
extent  is  assumed  to  be  Infinite. 

This  report  presents  a  method  for  satisfying  these  two  requirements  by 
coupling  two  solution  methods.  In  the  finite  element  method  used  to  provide 
a  detailed  analysis  of  the  inner  region,  the  inner  region  is  divided  into  a 
number  of  elements  so  that  local  changes  in  material  properties  and  local 
sources  can  be  represented  in  the  elements.  The  boundary  integral  equation 
method,  used  to  provide  the  analysis  of  the  outer  infinite  region,  reduces  the 
problem  of  solving  a  field  problem  in  the  external  region  to  a  much  smaller 
problem  on  just  the  Interface  between  the  inner  and  outer  regions.  The  prob¬ 
lem  is  discretized  by  dividing  the  Interface  into  a  number  of  segments.  A 
requirement  of  using  the  boundary  Integral  equation  method  in  the  form  de¬ 
scribed  in  this  report  is  that  the  outer  region  can  be  assumed  to  be  free  of 
sources  (acoustic,  magnetic,  etc.)  and  to  have  uniform  material  properties. 

This  report  avoids  the  complications  of  a  general  development  by  addres¬ 
sing  a  specific  problem,  namely  steady  state  heat  flow  in  two  dimensions.  Two 
formulations  of  the  coupling  are  given.  In  one  formulation  an  unsymmetric 
system  of  equations  derived  by  the  boundary  integral  equation  method  is  coupled 
to  the  symmetric  finite  element  system  of  equations,  producing  an  unsymmetric 
system  of  equations  to  be  solved.  In  the  other  formulation,  with  additional 
effort,  a  symmetric  system  is  derived  by  the  boundary  Integral  equation  method, 
which,  when  coupled  to  the  finite  element  equations,  produces  a  symmetric  sys¬ 
tem  of  equations  to  be  solved. 

A  sample  problem  was  solved  using  each  method.  In  this  sample  problem 
the  temperature  distribution  was  computed  for  steady  state  heat  conduction  in 
two-dimensions.  Both  methods  produced  good  results,  although  the  symmetric 
formulation  may  be  the  better. 
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II.  THE  PHYSICAL  PROBLEM  AND  ITS  MATHEMATICAL  FORMULATION 
In  considering  steady  state  heat  flow  in  two  dimensions.,  we  will 
determine  numerically  the  temperature  distribution  in  a  thin  slice  of  material 
of  infinite  extent  and  uniform  thermal  properties,  that  has  heat  sources  except 
in  a  bounded  region  I  (Figure  1).  Inside  this  inner  region  I,  the  thermal  con¬ 
ductivity  and  capacity  may  vary,  and  there  may  be  heat  sources  or  holes. 

Either  the  temperature  or  the  heat  flux  is  specified  on  the  boundary  of  the 
holes.  All  points  outside  the  inner  region  will  be  called  E,  the  outer  reg¬ 
ion. 


Figure  1  -  Diagram  of  the  Inner  and  Outer  Regions 

The  physical  problem  is  described  analytically  by  partial  differential 
equations  that  hold  in  the  two  regions  and  by  the  prescribed  data  on  boun¬ 
daries  of  the  regions.  The  Interface  between  the  inner  and  outer  regions  is 
called  F  and  the  boundaries  of  any  holes  in  the  inner  region  are  called  Fj. 

The  temperature  distribution  u  is  governed  by  the  differential  equations 

V^u  *  f,  in  I  (la) 

V^u  •  0,  in  E  (lb) 
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where  f  is  the  distribution  of  heat  sources  in  I.  On  the  boundaries  of  holes 
in  I,  the  temperature  and  heat  fluxes  satisfy  the  boundary  conditions 

u  =  g,  on  F  (2a) 


IS  =  h’  on  Fh  <2b> 

The  boundary  of  the  holes  is  composed  of  the  disjoint  parts  F  and  F.  on  which 

8  h 

are  prescribed  the  temperature  and  the  heat  flux,  respectively.  The  deriva¬ 
tive  in  the  direction  of  the  normal  vector  n  pointing  outward  from  the  region 

is  3/3n.  With  only  a  little  extra  complication,  we  could  also  consider  boun- 

3u 

dary  conditions  of  the  form  Au  +  B  7-  ■  C,  which  can  be  used  to  represent 

1*  dn 

convection  at  the  boundary. 


III.  THE  FINITE  ELEMENT  METHOD 

The  finite  element  method  is  used  to  solve  for  the  temperature  distri¬ 
bution  in  the  inner  region  which  has  the  more  complex  configuration  of  thermal 
properties.  To  apply  the  finite  element  method,  the  region  is  divided  into 
elements  that  are  small  enough  to  accommodate  changes  in  thermal  properties 
and  the  distribution  of  heat  sources  and  to  adequately  follow  the  shape  of  the 
boundaries.  In  practice  a  set  of  grid  points,  called  nodes,  is  defined  in  the 
region;  these  nodes  determine  the  sides  of  each  element.  Techniques  of  finite 

element  modeling  and  details  of  the  method  are  not  covered  in  this  report; 

2 

Zienklewicz  provides  a  comprehensive  treatment. 

The  finite  element  method  produces  a  matrix  equation  that  relates  the 
temperature  distribution  in  the  region  with  sources  of  heat  in  and  on  the 
boundary  of  the  region  and  properties  of  the  material  in  the  region.  This 
matrix  equation  is  of  the  form 

Ku  «  f  (3) 

where  K  is  a  matrix  determined  by  the  configuration  of  the  finite  elements  and 
the  thermal  properties  of  the  material,  u  is  the  vector  of  nodal  temperatures, 

*A  complete  listing  of  refer  ..  ea  is  given  on  page  29. 
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and  f  is  a  vector  of  nodal  heat  sources  determined  by  sources  inside  the  reg¬ 
ion  and  heat  fluxes  on  the  boundary. 

The  matrix  equation  is  modified  by  boundary  conditions.  For  a  node  at 
which  the  temperature  is  known,  the  component  of  u  for  the  node  is  constrained 
to  equal  the  known  value,  and  for  the  node  at  which  3u/3n  is  known,  terms  are 
added  to  the  source  vector  f.  On  the  interface  between  the  inner  and  outer 
regions,  the  boundary  values  serve  to  couple  the  solutions  for  the  two  regions. 
The  way  that  this  coupling  is  done  will  be  described  in  Section  V. 

IV.  THE  BOUNDARY  INTEGRAL  EQUATION  METHOD 
The  boundary  integral  equation  method  is  used  in  the  coupled  problem  to 
account  for  the  effects  of  the  infinite  outer  region  and  to  compute  the  tem¬ 
perature  at  selected  points  in  the  region.  The  method  reduces  the  problem  of 
solving  Equation  (lb)  for  all  points  in  the  region  to  the  problem  of  solving 
an  Integral  equation  on  the  boundary  of  the  region;  this  is  especially  help¬ 
ful  in  treating  the  infinite  outer  region.  Use  of  the  method  will  first  be 
described  for  a  bounded  region  and  will  then  be  extended  to  an  unbounded  region 
outside  a  curve.  This  description  is  offered  to  assist  in  understanding  the 

method  and  will  not  be  given  in  detail.  Improvements  and  refinements  to  the 

2  3  4 

method  have  been  described  by  Zienkiewicz,  Butterfield,  and  Brebbia. 

THE  INTERIOR  BOUNDARY  INTEGRAL  EQUATION  METHOD 

In  the  context  of  two-dimensional  steady  state  heat  flow,  the  problem  is 
to  determine  the  temperature  distribution  in  the  region  D  (see  Figure  2)  in¬ 
side  a  curve  F  when  the  temperature  is  known  on  part  of  the  boundary  F^  and 

the  heat  flow  is  known  on  the  remainder  of  the  boundary  F.  .  That  is,  the  dif- 

n 

ferential  equation.  Equation  (lb),  with  the  boundary  conditions.  Equations 
(2),  is  to  be  solved  by  reducing  the  problem  to  solving  an  integral  equation 
on  F. 

To  obtain  the  boundary  Integral  equation  we  begin  with  a  form  of  Green's 
theorem  for  the  region  D: 


in  which  V2  is  the  Laplace  operator  (3/3x)2  +  (3/3y)2  and  3/3n  is  the  deriva¬ 
tive  in  the  direction  of  the  normal  pointing  outward  from  D.  We  will  derive 
the  needed  integral  equations  by  substituting  two  particular  functions  for  the 
functions  v  and  w  in  Green's  theorem.  Equation  (lb)  and  the  boundary  condi¬ 
tions  (Equations  (2))  fora  a  well  posed  problem  having  a  continuous  solution  u 
which  we  substitute  for  w.  For  v,  we  substitute  the  function 
G(p,g)  ■  -  -|^-log  |p~q|  »  where  p  and  q  are  points  in  D;  p  is  the  point  at  which 
the  value  of  the  solution  u  is  obtained  and  q  is  the  variable  of  integration. 
This  form  of  G(p,q)  is  the  fundamental  solution  for  Equation  (1)  in  two- 
dimensions.  With  these  substitutions  Green's  theorem  becomes  (see  Appendix  A) 


u (p)  -  f( u(q)  ~  G(p,q)  )  ds(q)  (4) 

F 

For  any  point  p  in  D,  the  functions  G(p,q)  and  3G(p,q)/3n  are  explicity  known 
for  q  on  the  curve  F.  If  both  boundary  conditions  u  and  3u/3n  were  known  at 
all  points  on  the  boundary,  u(p)  could  be  computed  by  a  straightforward  inte¬ 
gration.  However,  for  a  well  posed  problem  exactly  one  of  the  two  boundary 
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conditions  must  be  specified  on  any  interval  of  the  boundary  so  in  the  context 
of  Equation  (4)  the  other  boundary  values  are  missing  on  this  interval.  The 
next  step  (which  is  the  first  step  in  the  numerical  implementation)  is  to 
obtain  the  missing  boundary  data. 

To  obtain  this  missing  data ,  we  move  the  point  p  to  the  boundary.  This 
must  be  done  carefully  (see  Appendix  B),  because  for  p  on  F,  the  integrand 
becomes  singular,  although  integrable.  For  p  on  F  Equation  (4)  becomes 

J  u(p)  =■  y^u(q)  3G^-^  ~  G(p,q)  )  ds(q)  (5) 

F 

This  is  the  boundary  integral  equation.  If  u  is  known  on  the  entire  boundary, 
that  Is,  if  we  have  a  Ln'richlet  problem,  3u/3n  is  unknown  and  the  boundary 
integral  equation  can  be  written 

/G(p,q)  |^ds(q)  =  /u(q)  ds(q)  -  \  u(p) 

F  F 

The  right  side  and  G(p,q)  are  explicitly  known.  This  integral  equation  is 
solved  for  3u/3n,  which  is  then  substituted  in  Equation  (4)  to  yield  tempera¬ 
tures  at  any  point  p.  In  a  similar  manner  if  3u/3n  is  known  on  the  boundary 
or  if  u  is  known  on  part  of  the  boundary  and  3u/3n  is  known  on  the  remainder, 
the  boundary  integral  equation  can  be  used  to  obtain  the  missing  boundary 
values . 

THE  EXTERIOR  BOUNDARY  INTEGRAL  METHOD 

To  solve  the  heat  conduction  problem  in  an  infinite  region,  we  will 
couple  the  finite  element  solution  in  a  bounded  part  of  the  region  to  the 
boundary  integral  solution  outside  of  the  bounded  part.  Therefore,  a  boundary 
Integral  solution  for  the  region  outside  of  a  curve  F  is  needed. 

To  compute  the  solution  u  to  Equation  (lb)  at  a  point  p  in  the  region 
outside  a  closed  curve,  we  begin  by  enclosing  both  the  point  p  and  the  curve  F 
in  a  circle  S  having  radius  R  (Figure  3).  Connecting  the  circle  and  the  curve 
F  with  a  cut  C  produces  a  region  to  which  the  preceding  work  applies.  The 
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integral  equation  for  points  inside  the  resulting  contour  and  for  points  on 
the  curve  F  takes  the  form 

au<p>  -  [/+/+/+/]  (G(p,q)  Mal.u(q)  ds(q) 

F  c+  C_  s 

where  a  *  1  if  p  is  inside  the  contour  and  a  *  1/2  if  p  is  on  the  curve  F  and 
F  has  a  smoothly  varying  tangent  at  p.  (Also  see  Appendix  A.) 


Figure  3  -  Contours  Used  to  Solve  the  Exterior  Problem 

In  physical  problems,  the  values  of  temperature,  sound  intensity,  magnetic 
field  strength,  etc.  go  to  zero  at  points  far  from  any  sources.  In  mathemati¬ 
cal  terms,  this  means  we  can  restrict  our  attention  to  functions  that  satisfy 
the  condition  u(q)  *  0(1/ q)  as  q  -*■  °°  .  Then  the  preceding  integral  con¬ 
verges  to  an  integral  over  the  curve  F  as  the  radius  R  of  the  circle  S  becomes 
infinite.  The  convergence  follows  from  two  facts:  (1)  The  integral  over  C+ 
equals  the  integral  over  C_  except  that  the  normal  vector  n  points  in  the 
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opposite  direction.  Thus  one  integrand  is  the  negative  of  the  other  and  the 

integrals  cancel.  (2)  Since  on  the  circle  S,  u(q)  *  0(1/R),  3u(q)/3n  « 

3u(q)/3r  *  0(1/R2)  and  lim  1°s  -  -  0, 

R-*»  R 

R i  log|p-qi  -  u(q)  log  I p-q  I  )  ds(q)  -  0 

S 

Thus  Equations  (4)  and  (5)  hold  for  points  p  outside  the  curve  F  if  the  noraal 
vector  n  points  inward,  that  is,  away  from  the  exterior  region. 


THE  COMPUTATIONAL  SOLUTION 

Equations  (4)  and  (5)  define  a  solution  to  the  boundary  integral  equation 
problem,  but  except  for  a  few  specially  shaped  boundaries  explicit  solutions 
are  not  available.  Therefore  a  means  of  obtaining  approximate  numerical  sol¬ 
utions  must  be  used.  To  avoid  unnecessary  complications  an  approach  in  which 
the  temperature  and  heat  flux  on  the  boundary  are  approximated  by  linear 
functions  is  used  here.  A  more  flexible  and  accurate  method,  which  uses  a  rep¬ 
resentation  of  the  boundary  segments  similar  to  that  used  in  isoparametric 
finite  elements,  is  described  by  Fairweather  et  al.* 

Equations  (4)  and  (5)  are  discretized  by  dividing  the  curve  F  into  k  seg¬ 
ments  Fj,  ...,  in  such  a  way  that  on  any  interval  F^  u  and  3u/3n  can  be 
approximated  by  linear  functions  and  that  on  any  one  Interval  either  u  or  3u/3n 

f 

is  given  as  boundary  data.  Let  u^  and  u^  be  the  approximations  to  the  values 
of  u  and  3u/3n  on  the  first  end  point  s^  of  the  segment  F^.  The  end  points 
are  ordered  by  traversing  the  boundary  in  the  direction  that  keeps  the  ex¬ 
terior  of  the  boundary  integral  region  on  the  left.  Let  N^(s)  be  the  linear 
function  on  the  segments  F^_j  and  F^  that  equals  unity  at  s^,  equals  zero  at 
8i+l  ^at  S2  and  8k*  *“1),  and  equal s  zero  on  all  other  segments  of  the 

boundary.  Then  the  linear  approximations  are 

k 

u(s)  «  E  u  N  (s)  (6a) 

i-i  1  1 


3u(s) 


k 

*  E 


u  N  (s) 


(6b) 


Substituting  these  approximations  into  Equations  (4)  and  (5)  produces  the 
discretized  equation 


au 


j 


k 

+  Z 
i-1 


3G(p^ ,q) 
3n 


ds(q) 


(7) 


where  a  -  1/2  if  p^  is  on  the  curve  F  and  F  has  a  continuously  turning  tangent 
at  Pj  (see  Appendix  A).  For  p^  in  the  interior  of  the  boundary  integral 
region,  a  ■  1, 

Since  G(p,q)  and  N^(s)  are  known  explicitly,  numerical  quadrature  gives 
numerical  values  for  the  matrix  components 

3G(p .q) 


JJ 

3n 


ij 


-  /  Ni(s) 

Fi+Fi-1 

/  Ni(s)  G(p,q)  ds(q) 


F  +F 
1  i-1 


For  i= j ,  the  integrand  is  singular,  but  its  value  can  be  found  by  an  analytic 
approximation.  The  matrix  A  is  defined  by  A„  =  A^  if  i  =  j,  and 

=  1/2  +  A^.  Combining  these  matrix  components  with  the  integral  Equation 
(7)  produces  the  matrix  equation 

Au  =  Bu 


where  u  =  (u. 


.j,  ...,  u^)  and  u  -  (Uj,  •••»  u^)  ate  the  vectors  of  temperatures 
and  normal  derivatives  at  the  end  points  of  the  segments.  This  is  the  matrix 
form  of  the  boundary  integral  equation. 

Recall  that  at  each  point  of  the  boundary  either  u  or  3u/3n  is  known  and 
the  other  is  to  be  determined,  so  that  exactly  half  of  the  components 


....  V  u1# 


,u^  are  known.  The  remainder  are  determined  by  exchanging 


r 

corresponding  columns  between  the  left  and  right  sides  of  the  matrix  equation 
as  necessary  to  move  the  unknown  components  to  the  left  side  and  the  known 
components  to  the  right.  This  procedure  yields  a  matrix  equation  of  the  form 


Mx  -  b 
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4. 


that  is  solved  for  the  unknown  boundary  values.  When  these  boundary  values 
have  been  computed,  all  the  boundary  values  are  known  at  all  the  segment  end 
points  on  the  boundary  F. 

Now  that  both  u  and  3u/3n  are  known  on  the  boundary,  we  can  use  Equation 
(7)  to  determine  the  temperature  at  any  point  p  that  is  not  on  the  boundary. 
This  time  we  use  the  matrix  components  A^  and  B^  to  get  the  equation 


u(Pj) 


k 

-  E 

i-1 


Vij 


k 

+  E 
i-1 


Ui  ®ij 


The  components  A^  and  have  the  same  form  for  p^  not  on  the  boundary  as 
the  corresponding  components  have  for  p^  on  the  boundary,  but  in  general,  they 
are  not  equal  to  any  of  these  components  and,  in  fact,  must  be  recomputed  for 
each  point  that  is  not  on  the  boundary. 


V.  COUPLING  THE  BOUNDARY  INTEGRAL  EQUATION 
AND  FINITE  ELEMENT  SOLUTIONS 

He  wish  to  solve  a  heat  flow  problem  using  a  finite  element  representation 
in  a  bounded  region  and  a  boundary  integral  equation  representation  of  the  sur¬ 
rounding  infinite  region.  The  individual  methods  have  been  described;  it 
remains  to  couple  the  solutions.  Two  formulations  of  the  coupling  will  be 
described.  The  first  formulation,  unsyrametrlc  coupling,  has  the  advantage  that 
the  boundary  integral  matrices  A  and  B  given  in  Equation  (8)  can  be  used  dir¬ 
ectly  without  the  need  of  any  matrix  inversion.  It  has  the  disadvantage  that 
the  coupling  introduces  unsymmetric  matrices  which  cause  the  entire  linear  sys¬ 
tem  to  require  an  unsymmetric  solution  method.  This  lack  of  symmetry  increases 
the  effort  needed  to  solve  the  system  of  equations.  The  advantage  of  the  sec¬ 
ond  formulation,  symmetric  coupling,  is  that  the  coupled  linear  system  remains 
symmetric,  so  faster  solution  routines  can  be  used  to  solve  the  system  of 
equations.  The  disadvantage  of  this  method  is  that  the  boundary  integral 
matrix  B  must  be  inverted.  If  this  matrix  is  small  enough  compared  to  the 
finite  element  matrix,  the  effort  saved  in  avoiding  an  unsymmetric  solution  may 
make  the  effort  of  inverting  the  matrix  worthwhile.  Further  details  and  some 
variants  of  these  coupling  formulations  are  described  by  Kelley  et  al.^ 
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UNSYMMETRIC  COUPLING 

The  unsymmetrlc  coupling  formulation  can  be  thought  of  as  a  finite  element 
solution  for  the  inner  region  with  special  boundary  conditions  applied  at  the 
interface.  The  boundary  integral  equation  method  produces  a  relation  between 
the  temperature  and  heat  flux  on  the  interface  that  provides  the  required 
boundary  conditions. 

The  finite  element  matrix  Equation  (3)  is  written  expressing  the  source 
term  as  the  irn  of  two  terms.  The  first  term  represents  the  heat  flux  through 
the  Interface  Into  the  inner  region  due  to  the  temperature  distribution  In  the 
outer  region.  The  equation  is 

Ku  «  f  +  f 
c 

The  i-th  component  of  fc  is 


(tc>i  -  / »  1  ds 

F 

where  the  N  are  the  linear  functions  (shape  functions)  defined  on  page  9 
*  2 

(see  Zienkiewica  ).  The  derivative  9u/9n  is  also  expanded  in  terms  of  the 
(Equation  (6b))  and  the  preceding  equation  becomes 


(tc>l 


k 

E 

i-1 


ds 


We  define  the  matrix  components 

cu  ’  /\Nj ds 
F 

and  write  the  finite  element  matrix  equation  in  the  form 

Ku  -  Cu'  +  f 


(9) 
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Now  we  need  to  relate  the  nodal  temperatures  u  with  the  nodal  fluxes  u'  in  this 
equation. 

The  boundary  integral  matrix  Equation  (8)  provides  the  needed  relation 
between  the  nodal  temperatures  and  fluxes  on  the  interface: 

.  _  I 

Au  ■  Bu 


The  preceding  two  matrix  equations  are  to  be  merged  to  form  the  coupled 

system.  Let  Uj . u^  and  uj,  ....  u^  be  the  nodal  temperatures  and  fluxes 

on  the  Interface,  and  let  u^j*  •  ••»  um  be  the  temperatures  at  the  remainder 
of  the  finite  element  nodes.  Merging  the  two  matrix  equations  produces  the 
following  coupled  finite  element-boundary  integral  system  of  equations: 


r  "1 

r  “I 

u 

l 

1 

i 

• 

i  c 

• 

1 

• 

1  1 

u 

1 

1 

k 

K  ! 

1 

u 

1 

k+1 

• 

!  o 

• 

1 

I 

• 

1 

u 

1 

\ 

m 

i  i 

t 

■  i 

u. 

i  i 

1 

■  i 

• 

A  I  0  J  -B 

• 

1  1 

• 

1  1 

"k  J 

(10) 


This  is  the  linear  system  of  equations  that  is  assembled  and  solved  by  a  com¬ 
puter  program  to  yield  temperatures  at  nodes  in  the  inner  region  and  tempera¬ 
tures  and  fluxes  on  the  Interface. 


SYMMETRIC  COUPLING 

The  symmetric  coupling  formulation  can  be  thought  of  as  a  more  unified 
approach.  It  is  derived  by  minimi  ring  an  energy  integral  for  both  the  inner 
and  outer  regions.  Both  the  finite  element  and  boundary  integral  equation 
methods  contribute  terms  to  the  energy  integral. 
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The  Integral  for  the  energy  In  the  Inner  and  outer  regions  and  for  energy 
arising  from  sources  on  the  boundaries  of  holes  in  the  inner  region  is 

J(u)  m  \  f  (?2u)dA  (?u>2dA  -  f  uhds 

E  1  \ 
m  1  J  *11  ds  +  \  J  CVu)2<iA  -  J  uhds 

F  IF. 

n 

2 

The  second  line  in  this  equation  follows  from  Green's  theorem  and  from  V  u  *  0 
in  the  region  E. 

To  obtain  a  matrix  form  of  the  energy  Integral  we  will  make  several  sub¬ 
stitutions  of  approximate  discretized  expressions.  The  second  term  in  the 
energy  integral  is  approximated  by  the  expression 


-T  - 
u  Ku 

2 

where  K  is  the  finite  element  stiffness  matrix  for  the  inner  region.  In  the 
first  and  third  terms,  we  substitute  the  expansions  (Equations  (6a),  (6b))  of 
u  and  3u/9n  in  terms  of  shape  functions  and  nodal  temperatures  and  fluxes, 
now  written  In  vector  notation  with  N  *  (N^,  ...,  N^) 

u  <55  Nu 

«  Nu' 

9n 

The  vector  of  fluxes  u'  can  be  expressed  in  terms  of  the  temperatures  u,  using 
the  boundary  Integral  matrix  Equation  (8)  solved  for  u* : 

u'  *  B  *Au  (H) 


so  that 


3u  _  -I  - 
Au 
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Also,  we  define  the  discretlr«d  source  vector 

»  -  (/  v»>») 

end  use  the  matrix  C  defined  by  Equation  (9).  With  these  approximations,  the 
energy  Integral  can  be  written  in  matrix  form 

-  IT  -1  -  1  -T  -  - 

J(u)  -  y  u  CB  Au  +  j  u  Ku  -  Hu 

Now  minimizing  the  energy  integral  J(u)  by  setting  its  variations  with  respect 
to  u^  equal  to  zero  yields  the  matrix  equation  for  the  temperature  vector  u 

|  (CTB"lA  +  (B_1A)TC)u  +  Ku  -  H  (12) 


5.  • 

I 

<4 


’i 


rT 

L  * 


This  is  the  symmetric  system  of  linear  equations  that  is  assembled  and  solved 
by  the  computer  program  to  yield  temperatures  at  nodes  in  the  inner  region  and 
on  the  Interface.  If  the  fluxes  are  needed  on  the  Interface,  they  are  easily 
obtained  by  saving  the  previously  computed  matrix  B  *A  and  usliq;  Equation  (11). 

The  problem  of  computing  the  coupled  finite  element-boundary  Integral 
solution  is  reduced  to  solving  either  of  the  systems  of  Equations  (10)  or  (12). 
Next,  we  develop  the  computer  methods  used  to  solve  these  systems. 

VI.  THE  COMPUTER  IMPLEMENTATION 

The  solution  method  described  in  the  preceding  section  is  Implemented 
using  two  existing  computer  programs.  A  boundary  Integral  equation  program 
which  employs  linear  approximations  is  used  to  form  the  matrices  A  and  B.^ 

This  program  is  easily  modified  to  save  the  matrices  in  a  form  that  can  be 
passed  to  the  finite  element  program.  The  computation  of  the  matrix  C 
(Equation  (9))  is  added  to  the  boundary  Integral  program.  The  NASTRAN  com¬ 
puter  program  solves  rather  general  finite  element  problems  and  therefore  can 
represent  the  inner  region  which  requires  detailed  representation.  In  addi¬ 
tion,  NASTRAN  will  easily  accept  and  manipulate  the  boundary  Integral  matrices 
and  incorporate  the  results  into  the  solution  of  the  finite  element  system  of 
linear  equations.  At  the  completion  of  the  finite  element  solution,  the 
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vectors  of  nodal  temperatures  and  heat  fluxes  are  passed  back  to  the  boundary 
Integral  program  for  use  in  computing  temperatures  required  in  the  outer 
region  using  Equation  (4). 

NASTRAN  contains  sequences  of  computations  which  are  executed  by  the 
program  to  carry  out  finite  element  solutions.  These  sequences,  called  rigid 
formats,  are  easily  altered  to  allow  the  introduction  of  variations  in  the 
finite  element  solution  procedures.  A  sequence  of  these  rigid  format  ALTERS  Is 
used  for  each  of  the  two  formulations  that  we  are  considering.  The  data  used 
by  NASTRAN  consist  of  the  finite  element  representation  of  the  inner  region, 
the  matrices  A,  B,  and  C,  and  the  list  of  rigid  format  ALTERS  which  give  the 
changes  used  to  incorporate  the  boundary  integral  matrices  into  the  finite  ele¬ 
ment  solution.  For  details  of  the  use  and  meaning  of  the  rigid  format  ALTERS 
see  the  NASTRAN  User's  Manual. ^ 

We  now  give  a  listing  and  description  of  the  rigid  format  ALTERS  used  for 
the  unsymmetrlc  coupling  formulation. 

ALTER  37 

INPUTT 2  /A,B,C,,/0/ll  $ 

PARTN  HKGG.PVL, /KSML , ,  $ 

MERGE,  ,,A,,PVS,/A0/1  $ 

MERGE,  ,C,,,,PVS/C0/1  $ 

MERGE  KSML,AO,CO,B,PVL,PVL/KCPL/1  $ 

EQUIV  RCPL.HKGG  $ 

ALTER  79,80 
DECOMP  HKLL/L.U  $ 

ALTER  93,99 
FBS  L.U.HPL/HULV  $ 

ALTER  108 

PARTN  HUGV, ,PVU/ ,PHI , ,/0  $ 

PARTN  HUGV, ,PVL/ ,DPH1, ,/0  $ 

OUT PUT 2  PHI,DPHI,,,//-1/11  $ 

ENDALTER 

The  INPUTT2  instruction  reads  the  matrices  A,  B,  and  C.  To  allow  room  for 
the  matrix  in  Equation  (10),  as  many  extra  degrees  of  freedom  are  generated 
as  there  are  nodes  on  the  interface,  so  the  finite  element  matrix  for  the  inner 


region  is  smaller  than  the  matrix  in  NASTRAN .  The  first  PARTN  instruction 
partitions  the  larger  matrix  HKGG  and  extracts  the  matrix  for  the  inner  region 
which  is  called  KSML.  The  three  MERGE  instructions  assemble  the  matrix  in 
Equation  (10)  and  call  it  KCPL.  The  EQUIV  instruction  equates  the  newly  gener¬ 
ated,  coupled  matrix  with  HKGG,  the  matrix  whose  equation  is  to  be  solved  as 
the  finite  element  computations  continue.  The  DECOMP  and  FBS  instructions 
ensure  that  unsymmetric  decomposition  routines  are  used  to  solve  the  system  of 
linear  equations.  After  the  finite  element  solution  is  complete,  the  two 
PARTN  and  the  0UTPUT2  Instructions  extract  the  temperatures  and  heat  fluxes  on 
the  Interface  and  save  them  for  use  by  the  boundary  Integral  equation  program. 

For  the  symmetric  coupling  formulation,  the  following  rigid  format  ALTERS 
are  used: 


ALTER  37 

INPUTT 2  /A,B,C,,/0/ll  $ 

SOLVE  B.A/E  $ 

MPYAD  C,E,/CTE/C,N,1  $ 

TRNSP  CTE/ETC  $ 

ADD  CTE,ETC/KPR/C,N,(.5,0.)/C,N,(.5,0.)  $ 

MERGE,  ,,,KPR,PVU,/KBI///C,N,6  $ 

ADD  HKGG.KBI/KCPL  $ 

EQUIV  KCPL, HKGG  $ 

ALTER  108 

PARTN  HUGV , , PVU/ , PHI , , /0  $ 

MPYAD  E,PHI,/DPHI/C,N,0  $ 

OUT PUT 2  PHI ,DPHI ,,,//— 1//11  $ 

ENDALTER 

The  INPUTT2  instruction  reads  the  matrices  A,  B,  and  C.  The  SOLVE,  MPYAD, 
and  TRNSP  Instructions  form  the  matrix  (CTB  *A  +  (B  *A)  C)/2  for  Equation  (12); 
this  matrix  is  called  KPR.  The  MERGE  Instruction  increases  the  dimension  of 
KPR  so  it  can  be  added  to  the  finite  element  matrix  HKGG  by  the  ADD  instruction 
which  forms  the  coupled  boundary  integral-finite  element  matrix  KCPL.  The 
EQUIV  instruction  sets  the  finite  element  matrix  HKGG  equal  to  the  coupled 
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matrix  KCrL  and  the  finite  element  solution  continues.  When  the  solution  is  com¬ 
plete,  the  PARTN  instruction  extracts  the  nodal  temperatures  on  the  interface, 
and  the  MPYAD  Instruction  produces  the  nodal  fluxes  using  the  matrix  E  computed 
earlier  by  the  SOLVE  instruction  according  to  Equation  (11).  The  0UTPUT2 
instruction  saves  the  nodal  temperatures  and  fluxes  for  use  by  the  boundary 
integral  program. 

If  temperatures  at  points  in  the  outer  region  are  desired,  the  boundary 
Integral  program  uses  the  nodal  temperatures  and  fluxes  on  the  interface  that 
are  saved  by  NASTRAN  to  compute  the  temperatures  according  to  Equation  (4). 

VII.  THE  SAMPLE  PROBLEM 

A  steady  state,  two-dimensional  heat  conduction  problem  was  solved  to 
Illustrate  the  coupled  finite  element  boundary  integral  method.  The  inner  and 
outer  regions  are  shown  in  Figure  4.  The  irregular  shape  was  chosen  to  avoid 
having  geometrical  symmetry  produce  better  results  than  would  otherwise  have 
been  obtained. 


Figure  4  -  Regions  and  Boundaries  for  the  Sample  Problem 
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The  sample  problem  Is  an  exterior  Dirichlet  problem.  The  temperature  is 
prescribed  on  the  Interior  boundary  F,  consistent  with  the  function 
u(r,0)  ■  cos  0/r,  where  r  is  the  distance  from  the  point  0  inside  the  boundary 
F,  and  0  is  the  polar  angle.  Since  u  is  a  known  solution  of  Equation  (1).  tem¬ 
perature  and  flux  values  can  be  calculated  at  any  point  to  compare  with  the 
computed  solution. 

Two  representations  of  the  inner  region  and  the  interface  in  terms  of 
finite  elements  and  boundary  segments  were  used;  these  representations  are 
shown  in  Figures  5a  and  5b.  The  grid  in  Figure  5b  has  twice  as  many  boundary 
segments  as  the  grid  in  Figure  5a,  but  except  for  adjustments  needed  along 
the  interface  to  accommodate  the  extra  points,  the  finite  element  grids  are 
the  same.  It  is  apparent  from  Figure  5  that  In  the  grid  with  fewer  boundary 
segments  the  length  of  the  segments  is  more  compatible  with  the  finite  element 
mesh  in  that  an  element  refinement  along  the  Interface  that  accommodates  only 
the  boundary  segments  is  not  required. 


Figure  5a  Figure  5b 

Figure  5  -  Boundary  Integral  and  Finite  Element  Grid  for  the  Sample  Problem 
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Results  from  the  sample  problem  are  tabulated  in  Appendix  B.  Four  trends 
can  be  seen  In  these  results: 

'  •  The  division  of  the  interface  into  segments  shorter  than  those 

consistent  with  the  finite  element  mesh  did  not  produce  more 
accurate  results  except  in  the  outer  region.  Even  there  the 
errors  are  not  significantly  smaller. 

'•  There  is  no  significant  difference  between  the  errors  from  the 
symmetric  and  unsymmetric  formulations,  except  in  the  outer 
region  where  the  errors  from  the  symmetric  formulation  average 
a  little  more  than  one-half  of  those  from  the  unsymmetric 
1  formulation. 

•  The  errors  in  both  the  temperatures  and  fluxes  on  the  interface 
become  larger  on  the  part  of  the  interface  where  there  are  fewer 
finite  elements  between  the  inner  boundary  and  the  interface. 

The  largest  errors  encountered  are  in  the  temperatures  and  fluxes 
on  the  interface.  Even  though  the  values  on  the  interface  are  sub¬ 
sequently  used  to  compute  temperatures  in  the  outer  region,  the 
errors  in  the  outer  region  are  quite  small. 

The  reasons  for  these  trends  have  not  been  studied  yet,  but  some  conjec¬ 
tures  can  be  made.  It  is  likely  that,  with  consistent  boundary  segment  lengths 
and  finite  element  meshes,  the  two  methods  produce  the  same  order  of  error,  so 
that  reducing  the  error  produced  by  one  part  of  the  procedure  does  not  improve 
the  overall  accuracy  of  the  coupled  method.  The  approximations  made  in  both 
the  boundary  integral  equation  and  finite  element  methods  are  the  same  for  the 
symmetric  and  unsymmetric  formulations,  so  it  would  be  expected  that  the 
results  would  he  nearly  equal.  Temperatures  in  the  outer  region  are  computed 
by  integrating  temperatures  and  fluxes  on  the  interface.  The  smoothing  char¬ 
acteristic  of  integration  tends  to  average  positive  and  negative  errors  on  the 
interface,  producing  results  in  the  outer  region  with  smaller  errors. 

The  results  of  this  admittedly  small  sample  indicate  that  the  method  pro¬ 
duces  generally  good  values  and  very  good  values  in  the  outer  region.  The  use 
of  boundary  segments  having  lengths  compatible  with  the  finite  element  mesh 
and  the  use  of  the  symmetric  formulation  are  recommended  when  the  coupled 
method  is  used. 
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APPENDIX  A  -  THE  BOUNDARY  INTEGRAL  EOUATION 
For  u  the  solution  of  the  differential  equation.  Equation  (lb),  that 
satisfies  the  boundary  conditions,  Equation  (2),  we  derive  Equation  (4)  and 
Equation  (5)  for  p  in  the  region  D  and  for  p  on  the  boundary  of  D,  respec¬ 
tively.  These  derivations  are  well  known  but  not  often  explicitly  written. 
Since  the  location  of  the  interface  between  the  inner  and  outer  regions  can  be 
chosen  so  that  u  is  harmonic  in  a  domain  that  extends  a  finite  distance  beyond 
D,  u  can  be  assumed  to  be  continuous  and  du/dn  to  be  bounded  in  and  on  the 

g 

boundary  of  D. 

For  p  in  D,  log  |p-q|  as  a  function  of  q  satisfies  the  Laplace  equation 

except  at  the  point  p.  A  small  disk,  with  center  p  and  radius  r  bounded  by 

the  circle  is  removed  from  D,  leaving  the  modified  region  Df  (Figure  6). 

On  D  ,  Green's  theorem  holds 
r 


(13) 


Figure  6  -  The  Modified  Region  Df  for  an  Interior  Point 

In  the  modified  region  Dr,  both  u  and  log  | p— q |  are  harmonic,  so  the  right 
side  in  the  preceding  equation  is  zero.  In  a  local  coordinate  system  with 
origin  at  p,  the  integral  over  K,  becomes 

W  (  -  Ir  108  r  +  r  )  r  d9 

0 
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! 

since  on  Kf,  jp— q |  -  r  and 
'  In  the  limit  as  r  0, 


and  Equation  (13)  becomes  Equation  (4): 

u(p)  -/  (U  3^  lo8 1  P-q  I  -  log|p-q|)  ds(q) 

F 

If  p  is  on  the  boundary  of  D  and  the  tangents  at  p  form  the  angle  a 
(Figure  7),  a  similar  derivation  provides  a  slightly  different  equation. 

This  time  we  remove  a  small  sector  with  vertex  p  and  radius  r  and  bounded  by 
the  arc  Af,  again  leaving  a  modified  region  Dr.  In  the  same  manner  as  before, 
taking  a  limit  and  applying  Green's  theorem  yields 

a 

lim  1 /*  /  3u  ,  u\  a.. 

r-0  2 VJ  ("arlo8  r  +7)r  d0  u(p) 

0 

and 

77  U<P>  m  J  log  |  p-q  |  +  u  log  |  p-q  |  j  ds(q) 

F 

This  result  is  more  general  than  Equation  (5),  which  is  obtained  for  the  case 
of  a  continuously  turning  tangent  at  p,  that  is,  a  -  it  . 


3/3n  «  -9/3r. 

2* 

hr!  77  r  log  r  d9  0 
0  2* 

h  f  u  de  -  u(p> 


Figure  7  -  The  Modified  Region  D  for  a  Boundary  Point 


APPENDIX  B  -  TABULATED  RESULTS  FROM  THE  SAMPLE  PROBLEM 
Tables  1  through  6  list  temperatures,  fluxes,  and  errors  computed  for  the 
sample  problem.  These  values  were  computed  using  the  symmetric  and  unsym- 
metric  formulations,  and  for  each  formulation  using  meshes  with  20  and  40  grid 
points  on  the  interface  between  the  inner  and  outer  regions  (see  Figure  5). 
Table  1  gives  temperatures  for  a  few  points  lying  in  the  outer  region  (grid 
points  numbered  1  to  6)  and  in  the  inner  region  (grid  points  numbered  7  to 
20).  Table  3  gives  temperatures  for  points  lying  on  the  interface  (grid  points 
numbered  21  to  40).  The  normal  derivatives  or  heat  fluxes  on  the  interface  are 
given  in  Table  5.  Tables  2,  4,  and  6,  which  correspond  to  Tables  1,  3,  and  5 
contain  relative  errors  in  the  computed  temperatures  and  fluxes  as  percents  of 
the  maximum  exact  value  for  each  group.  Figure  8  shows  the  location  of  the 
grid  points  for  which  results  are  tabulated. 

The  following  acronyms  are  used  to  designate  the  columns  in  the  tables: 

GP  -  grid  point  number 
X,Y  -  grid  point  coordinates 
U  -  exact  temperatures 
DU  -  exact  fluxes 

US  -  temperature  from  symmetric  formulation 
UU  -  temperature  from  unsymmetrlc  formulation 
DUS  -  flux  from  symmetric  formulation 
DUU  -  flux  from  unsymmetrlc  formulation 
ES  -  temperature  error  from  symmetric  formulation 
EU  -  temperature  error  from  unsymmetrlc  formulation 
MS  -  flux  error  from  symmetric  formulation 
DEU  -  flux  error  from  unsymmetrlc  formulation 
20,40  -  refer  to  grids  with  20  or  40  boundary  segments 
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TABLE  1  -  TEMPERATURES  IN  THE  INNER  AND  OUTER  REGIONS 


GP 

X 

Y 

U 

US  40 

US  20 

UU40 

UU20 

1 

8.00 

0.0700 

0.0735 

0.0778 

0.0805 

2 

2.00 

0.2500 

0.2638 

0.2651 

0.2673 

0.2679 

3 

6.00 

6.00 

0.0833 

0.0939 

0.0971 

0.1011 

0.1035 

4 

-6.00 

6.00 

-0.0833 

-0.0805 

-0.0771 

-0.0730 

-0.0699 

5 

-6.00 

-6.00 

-0.0833 

-0.0809 

-0.0776 

-0.0737 

-0.0708 

6 

6.00 

-6.00 

0.0833 

0.0945 

0.0977 

0.1014 

0.1036 

7 

2.26 

-0.74 

0.3996 

0.4192 

0.4177 

0.4226 

0.4201 

8 

1.19 

0.7608 

0.7780 

0.7771 

0.7799 

0.7786 

9 

1.05 

1.45 

0.3284 

0.3459 

0.3440 

0.3482 

0.3460 

10 

0.73 

1.01 

0.4702 

0.4806 

0.4794 

0.4823 

0.4808 

11 

0.48 

0.66 

0.7168 

0.7120 

0.7114 

0.7130 

0.7122 

12 

-1.01 

0.73 

-0.6472 

-0.6715 

-0.6662 

-0.6699 

-0.6633 

13 

-0.66 

0.48 

-0.9866 

-0.9965 

-0.9942 

-0.9955 

-0.9927 

14 

-1.45 

-0.4520 

-0.4669 

-0.4723 

-0.4628 

-0.4666 

15 

-0.66 

-0.9866 

-1.0053 

-1.0049 

-1.0038 

-1.0028 

16 

-0.74 

-2.26 

-0.1298 

-0.1319 

-0.1322 

-0.1282 

-0.1284 

17 

-0.39 

-1.19 

-0.2472 

-0.2565 

-0.2570 

-0.2541 

-0.2542 

18 

1.89 

-2.61 

0.1825 

0.1945 

0.1953 

0.1985 

0.1985 

19 

1.40 

-1.93 

0.2470 

0.2592 

0.2594 

0.2626 

0.2622 

20 

0.73 

-1.01 

0.4702 

0.4864 

0.4861 

0.4884 

0.4879 

TABLE  2  -  TEMPERATURE  ERRORS  IN  THE  INNER  AND  OUTER  REGIONS  - 
PERCENT  OP  MAXIMUM  TEMPERATURE 


GP 

X 

Y 

ES40 

ES20 

EU40 

EU20 

n 

6.00 

8.00 

1.8 

iff 

2.00 

2.00 

■ 

1.8 

mm 

6.00 

6.00 

1.8 

if 

-6.00 

6.00 

0.6 

1.0 

1.4 

mm 

-6.00 

-6.00 

■ 

0.6 

1.0 

1.3 

Kfl 

6.00 

-6.00 

i.i 

1.5 

1.8 

2.1 

7 

-0.74 

2.0 

1.8 

2.3 

8 

-0.39 

1.7 

1.6 

1.9 

9 

1.05 

1.45 

1.8 

1.6 

2.0 

1.8 

mm 

0.73 

1.01 

1.1 

0.9 

1.2 

1.1 

mm 

0.48 

wnzm 

0.5 

0.5 

0.4 

0.5 

El 

-1.01 

■I 

2.5 

1.9 

2.3 

1.6 

Kg 

-0.66 

1.0 

0.8 

0.9 

SKKKB'iW  B 

-1.45 

-1.05 

1.5 

2.1 

1.1 

1.5 

11a 

-0.66 

1.9 

1.9 

1.7 

1.6 

mm 

-0.74 

-2.26 

0.2 

0.2 

0.2 

0.1 

17 

-0.39 

-1.19 

0.9 

1.0 

0.7 

18 

1.89 

-2.61 

1.2 

1.3 

1.6 

1.6 

19 

1.40 

-1.93 

1.2 

1.3 

1.6 

1.5 

20 

0.73 

-1.01 

1.6 

1.6 

1.8 

1.8 

TABLE  3  -  TEMPERATURES  ON  THE  INTERPACE 


u 

US  40 

US  20 

UU40 

UU20 

0.1998 

0.2158 

0.2136 

0.2207 

0.2174 

0.1225 

0.1358 

0.1371 

0.1407 

0.1406 

0.0570 

0.0670 

0.0699 

0.0717 

0.0731 

0.0000 

0.0055 

0.0088 

0.0103 

0.0127 

-0.0669 

-0.0670 

-0.0611 

-0.0628 

-0.0577 

-0.1670 

-0.1731 

-0.1677 

-0.1683 

-0.1628 

-0.2996 

-0.3111 

-0.3044 

-0.3065 

-0.2990 

-0.4755 

-0.5102 

-0.5089 

-0.5033 

-0.4993 

-0.6024 

-0.6355 

-0.6351 

-0.6299 

-0.6266 

-0.6058 

-0.6329 

-0.6350 

-0.6303 

-0.6315 

-0.5219 

-0.5420 

-0.5433 

-0.5402 

-0.5393 

-0.3606 

-0.3714 

-0.3610 

-0.3712 

-0.3602 

-0.1626 

-0.1672 

-0.1607 

-0.1669 

-0.1598 

0.0000 

0.0021 

0.0077 

0.0046 

0.0115 

0.1298 

0.1394 

0.1416 

0.1422 

0.1443 

0.2470 

0.2620 

0.2652 

0.2648 

0.2675 

0.3210 

0.3372 

0.3410 

0.3401 

0.3434 

0.3522 

0.3727 

0.3753 

0.3762 

0.3774 

0.3333 

0.3531 

0.3526 

0.3572 

0.3554 

0.2749 

0.2936 

0.2945 

0.2982 

0.2973 

GP 

X 

21 

3.28 

22 

2.82 

23 

1.67 

24 

0.00 

25 

-1.43 

26 

-2.07 

27 

-2.18 

28 

-1.90 

29 

-1.66 

30 

-1.49 

31 

-1.25 

32 

-0.96 

33 

-0.59 

34 

0.00 

35 

0.74 

36 

1.40 

37 

2.04 

38 

2.57 

39 

3.00 

40 

3.29 

TABLE  4  -  TEMPERATURE  ERRORS  IN  THE  INTERFACE  - 
PERCENT  OF  MAXIMUM  TEMPERATURE 


-2.38 

-3.88 

-5.15 

-5.30 

-4.39 

-2.85 

-1.59 

-0.62 

0.00 

0.49 

0.91 

1.32 

1.81 

2.25 

2.26 
1.93 
1.48 
0.83 
0.00 

-1.07 
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TABLE  5  -  HEAT  FLUXES  ON  THE  INTERFACE 


X 

Y 

DU 

DUS40 

DUS  20 

DUU40 

DUU20 

-2.38 

0.0276 

0.0365 

0.0208 

0.0373 

0.0209 

-3.88 

0.0065 

0.0110 

0.0070 

0.0109 

0.0062 

1.67 

-5.15 

0.0055 

0.0073 

0.0054 

0.0064 

0.0036 

0.00 

-5.30 

0.0089 

0.0074 

0.0074 

0.0064 

0.0071 

-1.43 

-4.39 

0.0138 

0.0093 

0.0161 

0.0069 

0.0138 

-2.07 

-2.85 

0.0040 

-0.0032 

0.0055 

-0.0034 

0.0056 

-2.18 

-1.59 

-0.0263 

-0.0268 

-0.0013 

-0.0276 

-0.0044 

-1.90 

-0.62 

-0.1339 

-0.2070 

-0.1671 

-0.1976 

-0.1609 

-1.66 

0.00 

-0.3410 

-0.4075 

-0.3637 

-0.3946 

-0.3561 

-1.49 

0.49 

-0.3996 

-0.4440 

-0.4380 

-0.4440 

-0.4485 

-1.25 

0.91 

-0.3156 

-0.3446 

-0.3689 

-0.3508 

-0.3638 

-0.96 

1.32 

-0.1077 

-0.1294 

-0.0804 

-0.1445 

-0.0880 

-0.59 

1.81 

0.0594 

0.0364 

0.0630 

0.0220 

0.0568 

0.00 

2.25 

0.0659 

0.0529 

0.0757 

0.0480 

0.0792 

0.74 

2.26 

0.0590 

0.0581 

0.0597 

0.0562 

0.0569 

1.40 

1.93 

0.1131 

0.1206 

0.1243 

0.1286 

0.1229 

2.04 

1.48 

0.1440 

0.1484 

0.1504 

0.1462 

0.1493 

2.57 

0.83 

0.1371 

0.1543 

0.1499 

0.1543 

0.1478 

3.00 

0.00 

0.1031 

0.1163 

0.1021 

0.1174 

0.1016 

3.29 

-1.07 

0.0597 

0.0724 

0.0652 

0.0740 

0.0635 

TABLE  6  -  HEAT  FLUX  ERRORS  ON  THE  INTERFACE  - 
PERCENT  OF  MAXIMUM  FLUX 


GP 

X 

Y 

DES40 

DES20 

DEU40 

DEU20 

21 

3.28 

-2.38 

1.7 

2.4 

1.7 

22 

2.82 

-3.88 

0.1 

1.1 

0.1 

23 

1.67 

-5.15 

0.5 

0.0 

0.2 

0.5 

24 

0.00 

-5.30 

0.4 

0.4 

0.6 

0.5 

25 

-1.43 

-4.39 

1.1 

0.6 

1.7 

0.0 

26 

-2.07 

-2.85 

1.8 

0.4 

1.9 

0.4 

27 

-2.18 

-1.59 

0.1 

6.3 

0.3 

5.5 

28 

-1.90 

-0.62 

18.3 

8.3 

15.9 

6.8 

29 

-1.66 

0.00 

16.6 

5.7 

13.4 

3.8 

30 

-1.49 

0.49 

11.1 

9.6 

11.1 

12.2 

31 

-1.25 

0.91 

7.3 

13.3 

8.8 

12.1 

32 

-0.96 

1.32 

5.4 

6.8 

4.9 

33 

-0.59 

1.81 

5.8 

0.9 

0.7 

34 

0.00 

2.25 

3.3 

2.5 

3.3 

35 

0.74 

2.26 

MUMS 

0.2 

0.7 

0.5 

36 

1.40 

1.93 

1.9 

2.8 

3.9 

2.5 

37 

2.04 

1.48 

1.1 

1.6 

0.6 

1.3 

38 

2.57 

0.83 

3.2 

4.3 

2.7 

39 

3.00 

0.00 

0.3 

3.6 

0.4 

40 

3.29 

-1.07 

1.4 

3.6 

1.0 
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